Setup

<<ProPac`

?ProPac`GeoTools`*

?ProPac`Dynamics`*

?MakeODEs

x[t]/.t->0

PhaseTrajectory[f_,x_,T_,s1_,s2_]:=Module[{t,Eqns,vars,ICs,Solver,TT},
		vars=Map[ToExpression[ToString[#]<>"["<>ToString[t]<>"]"]&,x];
  	ICs=Map[ToExpression[ToString[#]<>"[0]"]&,x];
		Eqns=MakeODEs[x,f,t];
	  (*Solver[ODE_,vars,s_,TT_]:=
		NDSolve[Join[Eqns,{ICs[[1]]==0,ICs[[2]]==0}],vars,{t,0,TT}];*)
		Solver=NDSolve[Join[Eqns,{ICs[[1]]==s1,ICs[[2]]==s2}],vars,{t,0,T},MaxSteps->60000];
		ParametricPlot[vars/.Solver[[1]],{t,0,T},Compiled->False,DisplayFunction->Identity]
	]

PhasePortrait1[f_List,x_List,T_,R_?MatrixQ]:=Module[{s1,s2,R1,R2},
		R1=R[[1]];
		R2=R[[2]];
		Flatten[Table[PhaseTrajectory[f,x,15,s1,s2],{s1,R1[[1]],R1[[2]],R1[[3]]},{s2,R2[[1]],R2[[2]],R2[[3]]}]]
		]

?PhasePortrait

Example: Damped Pendulum

f={x2,-Sin[x1]-x2/2};x={x1,x2};

graphs2=PhasePortrait[f,x,8,{{-20,20,.5},{-3,3,3}}];

Show[graphs2,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{-10,10},{-3,3}},DisplayFunction->$DisplayFunction,Frame->True,FrameStyle->Black,ImageSize->{360,215},AspectRatio->Full,AxesOrigin->Automatic]

Show[%47,AxesStyle->Gray]

Show[%17,AxesStyle->Gray]

Show[%19]

Example: Soft Spring

U=x1^2/(1+x1^2);

F=Simplify[D[U,x1]];

f={x2,- x2/2-F};x={x1,x2};

graphs3=PhasePortrait[f,x,8,{{-20,20,.5},{-3,3,3}}];

Show[graphs3,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{-5,5},{-2,2}},DisplayFunction->$DisplayFunction,Frame->True,FrameStyle->Black,ImageSize->{360,215},AspectRatio->Full,AxesOrigin->Automatic]

Show[%42,AxesStyle->Gray]

Example: Crankshaft

U=x1^2+x1^4/4;

F=Simplify[D[U,x1]];

f={x2,(-D[2-Cos[2 x1],x1] x2^2/2-x2/2 -F)/(2-Cos[2 x1])};x={x1,x2};

graphs2=PhasePortrait[f,x,15,{{-6,6,2},{-5,5,5}}];

Show[graphs2,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{-4,4},{-4,4}},DisplayFunction->$DisplayFunction,Frame->True,FrameStyle->Black,ImageSize->{360,215},AspectRatio->Full,AxesOrigin->Automatic]

graphs3=Flatten[Table[PhaseTrajectory[f,x,15,s1,s2],{s1,-6,6,2},{s2,-5,5,5}]];

Show[graphs3,AxesLabel->{"x1","x2"},PlotRange->{{-4,4},{-4,4}},DisplayFunction->$DisplayFunction]

Hamiltonian

U=x1^2+x1^4/4;T=(2-Cos[2 x1]) x2^2/2;

sol=Solve[p==D[T,x2],x2]

H=Simplify[(T+U)/.sol[[1]]]

R=c x2^2/.sol[[1]]

Simplify[{D[H,p],-D[H,x1]-D[R,p]}]

f={D[H,p],-D[H,x1]-D[R,p]}/.c->1/2;x={x1,p}

graphs3=PhasePortrait[f,x,15,{{-6,6,2},{-5,5,5}}];

Show[graphs3,AxesLabel->{"x1","p"},PlotRange->{{-4,4},{-4,4}},DisplayFunction->$DisplayFunction]

Example: Van der Pol

f={x2,-x1+.8 (1-x1^2)x2};x={x1,x2};

graphs3=PhasePortrait[f,x,15,{{-6,6,.5},{-5,5,5}}];

Show[graphs3,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotRange->{{-4,4},{-4,4}},DisplayFunction->$DisplayFunction,Frame->True,FrameStyle->Black,ImageSize->{360,215},AspectRatio->Full]

Show[%18,ImageSize->{360,215},AspectRatio->Full]

Show[%15,Frame->True,FrameStyle->Black]

Options[Graphics]

{AlignmentPoint->Center,AspectRatio->Automatic,Axes->False,AxesLabel->None,,AxesOrigin->Automatic,AxesStyle->{},Background->None,BaselinePosition->Automatic,BaseStyle->{},ColorOutput->Automatic,ContentSelectable->Automatic,CoordinatesToolOptions->Automatic,DisplayFunction:>$DisplayFunction,Epilog->{},FormatType:>TraditionalForm,Frame->False,FrameLabel->None,FrameStyle->{},FrameTicks->Automatic,FrameTicksStyle->{},GridLines->None,GridLinesStyle->{},ImageMargins->0.`,ImagePadding->All,ImageSize->Automatic,ImageSizeRaw->Automatic,LabelStyle->{},Method->Automatic,PlotLabel->None,PlotRange->All,PlotRangeClipping->False,PlotRangePadding->Automatic,PlotRegion->Automatic,PreserveImageOptions->Automatic,Prolog->{},RotateLabel->True,Ticks->Automatic,TicksStyle->{}}

Example: Zero Dynamics

<<ProPac`

f={x1^2,-x2};x={x1,x2};

graphs3=PhasePortrait[f,x,8,{{-0.2,0.1,.05},{-0.5,0.5,0.5}}];

Show[graphs3,AxesLabel->{"x,(W^C)","y,(W^S)"},PlotRange->{{-0.2,0.2},{-0.25,0.25}},DisplayFunction->$DisplayFunction,Frame->True,FrameStyle->Black,ImageSize->{360,215},AspectRatio->Full,AxesOrigin->Automatic,AxesStyle->Directive[Black,Thick]]

graphs=PhasePortrait[f,x,8,{{-.2,.1,.05},{-.5,.5,.5}}];

Show[graphs,AxesLabel->{"x","y"},PlotRange->{{-.2,.2},{-.5,.5}},DisplayFunction->$DisplayFunction]